city_code<-'qd' #目标城市
pty<-11    #产品类型
minmum_city<-30  #城市最小样本量设定
minmum_dist<-30  #行政区最小样本量设定
# date_base<-24109  #基期设定

library(RMySQL)
con <- dbConnect(MySQL(),host="192.168.3.139",
                 dbname="niek_160906",user="niek",password="FbWf8AKC")  
dbSendQuery(con,'set names utf8')

#抽取数据
source('R/load_spatial_data.R')
# price_sql<-paste0("select city_code,ymid,ha_code,proptype,saleprice,salebldgarea,salecount from ha_price where city_code=\'",city_code,"\'","and proptype=",pty)
# price_sale<-dbGetQuery(con,price_sql)
# info<-dbGetQuery(con,paste0("select * from ha_info where city_code=\'",city_code,"\'"))
# bldg<-dbGetQuery(con,paste0("select * from ha_bldg where city_code=\'",city_code,"\'"))
# phase<-dbGetQuery(con,paste0("select * from ha_phase where city_code=\'",city_code,"\'"))

#ha_price与ha_phase数据处理
library(dplyr)
price_sale<-subset(tabln.vec$ha_price,proptype==pty)[,c('ymid','ha_code','saleprice','salebldgarea')]%>%na.omit
phase<-tabln.vec$ha_phase%>%
    filter(!is.na(buildyear))%>%
    group_by(ha_code)%>%
    summarise(buildyear=round(mean(buildyear)))%>%as.data.frame

#关联phase数据--建筑年代
pp_sale<-merge(price_sale,phase,by='ha_code')

#关联bldg数据
# library(nnet)
# bldg2<-cbind(bldg,class.ind(bldg$bldg_type))
# pib_sale<-merge(pi_sale,bldg2[,c(-1,-3,-4)],by='ha_code',all.x=TRUE)

#关联POI与info数据
ppi_sale<-merge(pp_sale,ha_info.sp,by='ha_code')

# ############################################################################
# #基期与基点设定
# date_base<-24133  #基期设定为2011年1月
# city_base<-'bj'   #城市设定为青岛
# data_base<-subset(ppi_sale,ymid==date_base&city_code==city_base)
# myvar_city<-names(data_base) %in% c('ymid','city_code','ha_code','proptype','salecount','name','ha_cl_code','ha_cl_name','dist_code','dist_name','x','y','volume_rate','greening_rate')
# data_base1<-data_base[!myvar_city]
# 
# if(nrow(data_base1)>=300){data_base2<-data_base1[1,]}else{stop("该基期样本量不足300，需重新设定")}
# 
# len_base<-length(names(data_base2))
# for(ii in 1:len_base){
#     name_base<-names(data_base2)[ii]
#     data_base2[,name_base]<-mean(data_base1[,name_base],na.rm=T)
# }
# 
# fit_base<-lm(log(saleprice)~.,data = data_base1)
# price.pre_base<-predict(fit_base,newdata=data_base2)[[1]]%>%exp
# ############################################################################

#城市级别模型建立
cat(city_code,'model building...','\n')
myvar_city<-names(ppi_sale) %in% c('city_code','ha_code','proptype','salecount','name','ha_cl_code','ha_cl_name','dist_code','dist_name','x','y','volume_rate','greening_rate')
sale<-ppi_sale[!myvar_city]

#城市最小样本量设定
tab.date<-table(ppi_sale$ymid)%>%as.data.frame%>%subset(Freq>=minmum_city)

date<-sort(unique(tab.date$Var1))%>%as.character%>%as.numeric
len<-length(date)
fit.vec_city<-vector(mode='list',len)
names(fit.vec_city)<-date
for (i in 1:len){
    sale.new<-subset(sale,ymid==date[i])
    fit.vec_city[[i]]<-lm(log(saleprice)~.,data = sale.new[,names(sale.new)!='ymid'])
}

#计算城市HPI（拉氏指数）
cat('compute HPI of',city_code,'\n')
# data_base<-sale[1,-1]
# data_base1<-subset(sale,ymid==date[1]) #基期设定为最早日期
# 
# len_base<-length(names(data_base))
# for(ii in 1:len_base){
#     name_base<-names(data_base)[ii]
#     data_base[,name_base]<-mean(data_base1[,name_base],na.rm=T)
# }

result<-as.data.frame(cbind(dt<-rep(0,len),price.pre<-rep(0,len),r2<-rep(0,len),FBPI<-rep(0,len)))
names(result)<-c('DT','PRE','AdjR2','FBPI')
for(j in 1:len){
    price.pre<-predict(fit.vec_city[[j]],newdata=data_base2)
    fit.sum<-summary(fit.vec_city[[j]])
    fit.squ<-fit.sum$adj.r.squared
    result[j,1]<-date[j]
    result[j,2]<-exp(price.pre[[1]])
    result[j,3]<-fit.squ
    result[j,4]<-result[j,2]/price.pre_base*100
}


#行政区级别模型建立
myvar_dist<-names(ppi_sale) %in% c('city_code','ha_code','proptype','salecount','name','ha_cl_code','ha_cl_name','dist_name','x','y','volume_rate','greening_rate')
ppi_sale1<-ppi_sale[!myvar_dist]
dist.code<-as.vector(unique(ppi_sale1$dist_code))
len_dist<-length(dist.code)
fit.vec_dist<-vector(mode='list',len_dist)
names(fit.vec_dist)<-dist.code
for(m in 1:len_dist){
    dist_cd<-dist.code[m]
    cat(dist_cd,'model building...','\n')
    sale_dist<-subset(ppi_sale1,dist_code==dist_cd)
    
    #行政区最小样本量设定
    tab.date_dist<-table(sale_dist$ymid)%>%as.data.frame%>%subset(Freq>=minmum_dist)
    tab.date_dist1<-tab.date_dist$Var1%>%as.character%>%as.numeric
    date_dist<-sort(unique(tab.date_dist1))
    len_date_dist<-length(date_dist)
    
    if(len_date_dist==0) next  #如果每月样本量小于30，跳出本次循环
    
    fit.vec_dist1<-vector(mode='list',len_date_dist)
    names(fit.vec_dist1)<-date_dist
    for (n in 1:len_date_dist){
        sale_dist.new<-subset(sale_dist,ymid==date_dist[n])
        fit.vec_dist1[[n]]<-lm(log(saleprice)~.,data = sale_dist.new[,names(sale_dist.new)!='ymid' & names(sale_dist.new)!='dist_code'])
    }
    fit.vec_dist[[m]]<-fit.vec_dist1
}

#计算行政区HPI（拉氏指数）
result.vec_dist<-vector(mode='list',len_dist)
names(result.vec_dist)<-dist.code
for(pp in 1:len_dist){
    dist_cd<-dist.code[pp]
    cat('compute HPI of',dist.code[pp],'\n')
    sale_dist<-subset(ppi_sale1,dist_code==dist_cd)
    
    #行政区最小样本量设定
    tab.date_dist<-table(sale_dist$ymid)%>%as.data.frame%>%subset(Freq>=minmum_dist)
    tab.date_dist1<-tab.date_dist$Var1%>%as.character%>%as.numeric
    date_dist<-sort(unique(tab.date_dist1))
    len_date_dist<-length(date_dist)
    
    if(len_date_dist==0) next  #如果每月样本量小于30，跳出本次循环
    
    # data_base_dist<-sale_dist[1,names(sale_dist)!='dist_code'&names(sale_dist)!='ymid']
    # data_base_dist1<-subset(sale_dist,ymid==date_dist[1]) #基期设定为最早日期
    # 
    # len_base_dist<-length(names(data_base_dist))
    # for(qq in 1:len_base_dist){
    #     name_base_dist<-names(data_base_dist)[qq]
    #     data_base_dist[,name_base_dist]<-mean(data_base_dist1[,name_base_dist],na.rm=T)
    # }
    
    result_dist<-as.data.frame(cbind(dt<-rep(0,len_date_dist),price.pre<-rep(0,len_date_dist),r2<-rep(0,len_date_dist),FBPI<-rep(0,len_date_dist)))
    names(result_dist)<-c('DT','PRE','AdjR2','FBPI')
    
    for(jj in 1:len_date_dist){
        fit_dist<-fit.vec_dist[[pp]][[jj]]
        price.pre_dist<-predict(fit_dist,newdata=data_base2)
        fit_dist.sum<-summary(fit_dist)
        fit_dist.squ<-fit_dist.sum$adj.r.squared
        result_dist[jj,1]<-date_dist[jj]
        result_dist[jj,2]<-exp(price.pre_dist[[1]])
        result_dist[jj,3]<-fit_dist.squ
        result_dist[jj,4]<-result_dist[jj,2]/price.pre_base*100
    }
    result.vec_dist[[pp]]<-result_dist
}

#城市与行政区最终FBPI、PRE、AdjR2
result_FBPI<-result[,c('DT','FBPI')]
names(result_FBPI)<-c('DT',city_code)

result_PRE<-result[,c('DT','PRE')]
names(result_PRE)<-c('DT',city_code)

result_AdjR2<-result[,c('DT','AdjR2')]
names(result_AdjR2)<-c('DT',city_code)

for(xx in 1:len_dist){
    dist_cd<-dist.code[xx]
    result.vec_dist.new<-result.vec_dist[names(result.vec_dist)==dist_cd][[1]]
    
    if(is.null(result.vec_dist.new)) next  #如果结果为NULL，跳出本次循环
    
    result.vec_dist.new_FBPI<-result.vec_dist.new[,c('DT','FBPI')]
    names(result.vec_dist.new_FBPI)<-c('DT',dist_cd)
    result_FBPI<-merge(result_FBPI,result.vec_dist.new_FBPI,by='DT',all=TRUE)
    
    result.vec_dist.new_PRE<-result.vec_dist.new[,c('DT','PRE')]
    names(result.vec_dist.new_PRE)<-c('DT',dist_cd)
    result_PRE<-merge(result_PRE,result.vec_dist.new_PRE,by='DT',all=TRUE)
    
    result.vec_dist.new_AdjR2<-result.vec_dist.new[,c('DT','AdjR2')]
    names(result.vec_dist.new_AdjR2)<-c('DT',dist_cd)    
    result_AdjR2<-merge(result_AdjR2,result.vec_dist.new_AdjR2,by='DT',all=TRUE)
}

#保存结果
# file_name<-paste0(city_code,'_sale',pty)
# file.save<-paste0("~/mnt/HPI/",file_name,"new.csv")
# write.csv(result,file.save)


